"""

Computational Economics
05: Loops
http://johnstachurski.net/lectures/loops.html

"""


import math, random


"""
    Problem 1:

        * Get n from user and compute n! with a for loop
        * Recall that n! = n * (n - 1) * ... * 2 * 1
"""
prompt = 'Input a number and we will factorialize it: '
n = int(raw_input(prompt))
n_fact = 1

for i in range(n):
    n_fact *= n - i

print '%s! = %s' % (n, n_fact)



"""
    Problem 2: (Redo an earlier exercise using a for loop)

        * Simulate a draw from Y ~ Bin(n, p) using uniform rvs.
              o num of successes in n indep trials with success prob p
              o Get n and p from the user
"""
prompt = 'Input the number of trials: '
n = int(raw_input(prompt))

prompt = 'What is the % chance of success (1-99): '
p = float('%.2f' % (float(raw_input(prompt)) / 100))
successes = 0

for i in range(n):
    if (random.uniform(0, 1) < p):
        successes += 1

print 'num of trial successes (n=%s, p=%s): %s' % (n, p, successes)



p3 = """
    Problem 3:

        * Compute an approximation to pi using Monte Carlo
              o Note that if U is uniform on the unit square, then the probability U is in a subset B is equal to the area of B
              o And that if U_1,...,U_n are IID copies of U, then the fraction in B converges to the probability of B as n gets large
              o Finally, recall that for a circle, area = pi * radius^2
"""
print p3
print 'MY SOLUTION:'
plots = 100000
radius = 1

in_circle = 0
for n in range(plots):
    x, y = random.uniform(0, radius), random.uniform(0, radius)
    c = math.sqrt(x**2 + y**2)
    if c < radius:
        in_circle += 1

monte_carlo_fraction = '%.4f' % ((in_circle / float(n)) * 4)

print 'pi monte carlo (plots=%s, in_circle=%s): %s' % (plots, in_circle,
    monte_carlo_fraction)


print 'BOOK SOLUTION:'
from random import uniform
from math import sqrt

n = 100000

count = 0
for i in range(n):
    U, V = uniform(0, 1), uniform(0, 1)
    d = sqrt((U - 0.5)**2 + (V - 0.5)**2)
    if d < 0.5:
        count += 1

area_estimate = count / float(n)

print 'area estimate: %s' % (area_estimate * 4)  # dividing by radius**2


"""
    Problem 4:

        * Write a program which prints one random outcome of following game
              o 10 flips of an unbiased coin
              o If 3 consecutive heads occur, pays one dollar
              o If not, pays nothing
"""
flips = 10
results = []
for n in range(flips):
    result = (random.uniform(0,1) > 0.5) and 'H' or 'T'
    results.append(result)

paid = ''.join(results).find('HHH') != -1

print 'results %s paid: %s' % (results, paid)



"""
    Problem 5:

        * Use Monte Carlo to estimate average payoff of previous game
"""

# repeated trials
num_trials = 10000

def flips_paid():
    flips = 10
    results = []
    for n in range(flips):
        result = (random.uniform(0,1) > 0.5) and 'H' or 'T'
        results.append(result)

    return ''.join(results).find('HHH') != -1

winners = sum([flips_paid() for n in range(num_trials)])
print 'results paid in %s trials: %s%%' % (num_trials,
    '%.1f' % (100 * float(winners) / num_trials))



print '%s: ok' % (__file__)
